# MEDIPS batch processing
library(MEDIPS)
library(BSgenome.Hsapiens.UCSC.hg19)

getwd()

#get list of files in directory
files <- list.files(pattern=".txt")

#####################################
# beginning of pre-processing loop
#####################################

#loop through files, process and write results to individual RData files 
for (i in files) {
	sample <- gsub("_sorted.txt.RData","",i)
	# road alignments
	MEDIPS.SET <- MEDIPS.readAlignedSequences(BSgenome = "BSgenome.Hsapiens.UCSC.hg19", file = i)
	# reads are extended to 250bp, MEDIPS manual uses 400bp
	# 250bp is based on observation of mean insert size ~180bp
	# bin_size (50) could also be adjusted, another commonly used value for this parameter is 500
	MEDIPS.SET <- MEDIPS.genomeVector(data = MEDIPS.SET, bin_size = 50, extend = 250)
	if (which(files == i) == 1) {
		# create the genome vector for the first set and store results in variables
		MEDIPS.SET <- MEDIPS.getPositions(data = MEDIPS.SET, pattern = "CG")
		pattern_chr <- MEDIPS.SET@pattern_chr
		pattern_pos <- MEDIPS.SET@pattern_pos
		seq_pattern <- MEDIPS.SET@seq_pattern
		number_pattern <- MEDIPS.SET@number_pattern
		# calculate the CpG coupling factors across the genome for the first set and store results in variables
		# manual suggests 700
		MEDIPS.SET <- MEDIPS.couplingVector(data = MEDIPS.SET, fragmentLength = 500, func = "count")
		extend <- MEDIPS.SET@extend
		fragmentLength <- MEDIPS.SET@fragmentLength
		genome_CF <- MEDIPS.SET@genome_CF
		distFunction <- MEDIPS.SET@distFunction
		distFile <- MEDIPS.SET@distFile
	} else {
		# for subsequent sets the genome vector information can be reused from initial analysis
		MEDIPS.SET@pattern_chr <- pattern_chr
		MEDIPS.SET@pattern_pos <- pattern_pos
		MEDIPS.SET@seq_pattern <- seq_pattern
		MEDIPS.SET@number_pattern <- number_pattern
		# for subsequent sets the coupling factor information can be reused from initial analysis
		MEDIPS.SET@extend <- extend
		MEDIPS.SET@fragmentLength <- fragmentLength
		MEDIPS.SET@genome_CF <- genome_CF
		MEDIPS.SET@distFunction <- distFunction
		MEDIPS.SET@distFile <- distFile
	}

	#####################################
	# saturation analysis
	#####################################
	sa.MEDIPS.SET <- MEDIPS.saturationAnalysis(data = MEDIPS.SET, bin_size = 50, extend = 250, no_iterations = 10, no_random_iterations = 1)
	png(file = paste(sample,"_saturationAnalysis.png",sep=""), width=1024, height=1024)
	MEDIPS.plotSaturation(sa.MEDIPS.SET)
	dev.off()
	
		# collating saturation analysis data for all samples
		if (which(files == i) == 1) {
			sa.ALL <- as.matrix(cbind(sa.MEDIPS.SET$numberReads, sa.MEDIPS.SET$maxEstCor[2], sa.MEDIPS.SET$maxTruCor[2]))
			colnames(sa.ALL) <- c("numberReads", "maxEstCor", "maxTruCor")
			rownames(sa.ALL)[which(files == i)] <- sample
		} else {
			sa.ALL <- rbind(sa.ALL, cbind(sa.MEDIPS.SET$numberReads, sa.MEDIPS.SET$maxEstCor[2], sa.MEDIPS.SET$maxTruCor[2]))
			rownames(sa.ALL)[which(files == i)] <- sample
		}
	
		# create calibration curve
		MEDIPS.SET <- MEDIPS.calibrationCurve(data = MEDIPS.SET)
		png(file = paste(sample,"_calibrationCurve.png", sep=""), width=1024, height=1024)
		MEDIPS.plotCalibrationPlot(data = MEDIPS.SET, linearFit = T, xrange=100, main = paste("Sample ", sample, sep=""))
		dev.off()
	
	#####################################
	# perform coverages analysis, 
	# i.e. coverage of CpG across genome 
	#####################################
	ca.MEDIPS.SET <- MEDIPS.coverageAnalysis(data = MEDIPS.SET, extend = 250, no_iterations = 10)
	png(file = paste(sample,"_coverageAnalysis.png",sep=""), width=1024, height=1024)
	MEDIPS.plotCoverage(ca.MEDIPS.SET)
	dev.off()
	
		# collating coverage analysis results
		if (which(files == i) == 1) {
			ca.ALL.CpGs.covered <- as.matrix(rbind(ca.MEDIPS.SET$coveredPos[2,]))
			colnames(ca.ALL.CpGs.covered) <- c("x1","x2","x3","x4","x5","x10")
			rownames(ca.ALL.CpGs.covered)[which(files == i)] <- sample
			ca.ALL.CpGs.fraction <- as.matrix(rbind(ca.MEDIPS.SET$coveredPos[3,]))
			colnames(ca.ALL.CpGs.fraction) <- c("x1","x2","x3","x4","x5","x10")
			rownames(ca.ALL.CpGs.fraction)[which(files == i)] <- sample
		} else {
			ca.ALL.CpGs.covered <- rbind(ca.ALL.CpGs.covered, ca.MEDIPS.SET$coveredPos[2,])
			rownames(ca.ALL.CpGs.covered)[which(files == i)] <- sample
			ca.ALL.CpGs.fraction <- rbind(ca.ALL.CpGs.fraction, ca.MEDIPS.SET$coveredPos[3,])
			rownames(ca.ALL.CpGs.fraction)[which(files == i)] <- sample
		}
	
	#####################################
	# CpG enrichment analysis
	#####################################
	er.MEDIPS.SET = MEDIPS.CpGenrich(data = MEDIPS.SET, extend = 200)
	
		# collating enrichment analysis data for all samples
		if (which(files == i) == 1) {
			er.ALL <- as.matrix(rbind(unlist(er.MEDIPS.SET)))
			rownames(er.ALL)[which(files == i)] <- sample
		} else {
			er.ALL <- rbind(er.ALL, unlist(er.MEDIPS.SET))
			rownames(er.ALL)[which(files == i)] <- sample
		}
	
	######################################
	# normalization, CpG-density dependent
	######################################
	MEDIPS.SET = MEDIPS.normalize(data = MEDIPS.SET)

	######################################
	# saving processed data and 
	# remove objects before next iteration
	######################################
	save(list = c("MEDIPS.SET", "er.MEDIPS.SET", "sa.MEDIPS.SET", "ca.MEDIPS.SET"), file = paste(sample, ".RData", sep="")))
	rm(list = c("MEDIPS.SET", "er.MEDIPS.SET", "sa.MEDIPS.SET", "ca.MEDIPS.SET"))
	
	gc()
}

###########################################
# end of pre-processing loop
###########################################

###########################################
# saving QC analysis results of all samples
###########################################
save(list = c("er.ALL","sa.ALL","ca.ALL.CpGs.covered","ca.ALL.CpGs.fraction"), file = "AnalysisSummaries.RData")





